// A : atmic mass
// Mt : target mass
// Mwimp : wimp mass
// sigmawn : Xsec for wimp and nucleon

double c = 3.0e8; // [m/s]
double ck = 3.0e5; // [km/s]
double rho0 = 0.3; //[GeV/c2/cm3] WIMP energy density
double Mwimp = 100; // [GeV/c2] WIMP mass
double Mn = 0.9315; // [GeV/c2] average nucleon mass
double v0 = 220;    // [km/s] Sun's velocity from Galactic center
double vesc = 544;  // [km/s] Galactic escape velocity
//double vesc = 600;  // [km/s] Galactic escape velocity
double vear = 232;  // [km/s] Average earth velocity from Galactic center
double NA = 6.0221415e23;  // Avogadro's number
double pi = TMath::Pi();

double c1 = 0.751;      //factor change for each months
double c2 = 0.561;      //factor change for each months

double sigma0(double A, double Mt, double Mw, double sigmawn); 
double Mr(double M1, double M2); // reduced mass
double F(double Q, double A);

double dRdE(double ER ,double A , double sigmawn){
  // argon 

  double Mt = A * 0.932; // GeV/c2
  double sigmaA = pow(A,2) * pow( Mr(Mwimp,Mt)/Mr(Mwimp,Mn),2 ) * sigmawn;
  double vmin = sqrt( Mt*ER / (2*Mr(Mwimp,Mt)) );
  double Q = sqrt(2*Mt*ER)/c; // keV s/m

  double k0 = pow(pi*pow(v0,2),3/2);
  double kf = k0 * ( TMath::Erf(vesc/v0) - 2/sqrt(pi)*vesc/v0*exp( -pow(vesc/v0,2) ) );
  //  double R0 = 2 *NA * rho0/sqrt(pi)/A/Mwimp * sigmaA * v0;
  double R0 = (361/Mwimp/Mt) * (sigmaA) * (rho0/0.3) * (v0/220);   // GeV , pb , GeV/cm3 , km/s  ## Nishimura ph.D ##

  double E0 = Mwimp*pow((v0/ck),2)/2.;         // GeV*(km/s)^2   E(100GeV)=mc^2
  double rm = 4*Mwimp*Mt/pow(Mwimp+Mt,2); 

  //  double Term1 = sqrt(pi)*v0/(4*vear) * ( TMath::Erf((vmin+vear)/v0) - TMath::Erf((vmin-vear)/v0) );
  double Term1 = c1 * exp(-c2 * ER / E0 / rm);  // year average
  double Term2 = exp( -pow(vesc/v0,2) );
 
  //  double drate =  R0/E0/rm * exp(-ER/E0/rm) *1e-6;  // vesc = inf , vear = 0

  //  double drate = k0/kf * R0/E0/rm * ( Term1 - Term2 ) *1e-6; //* pow(F(Q,A,Mt),2);    
  double drate = k0/kf * R0/E0/rm * ( Term1 - Term2 ) *1e-6; //  vesc = vesc , vear = vear 

    //  return  R0/E0/rm * exp(-ER/E0/rm) *1e-6;
    return drate;
}


double dRdE2(double *x ,double *par){
  // argon 
  double ER = x[0]*1e-6; // keV
  double A = par[0];  // atomic mass
  double Mt = A * 0.932; // GeV/c2
  double sigmawn = par[1]; // pb

  double sigmaA = pow(A,2) * pow( Mr(Mwimp,Mt)/Mr(Mwimp,Mn),2 ) * sigmawn;
  double vmin = sqrt( Mt*ER / (2*Mr(Mwimp,Mt)) );
  double Q = sqrt(2*Mt*ER)/c; // keV s/m

  double k0 = pow(pi*pow(v0,2),3/2);
  double kf = k0 * ( TMath::Erf(vesc/v0) - 2/sqrt(pi)*vesc/v0*exp( -pow(vesc/v0,2) ) );
  //  double R0 = 2 *NA * rho0/sqrt(pi)/A/Mwimp * sigmaA * v0;
  double R0 = (361/Mwimp/Mt) * (sigmaA) * (rho0/0.3) * (v0/220);   // GeV , pb , GeV/cm3 , km/s

  double E0 = Mwimp*pow((v0/ck),2)/2.;         // GeV*(km/s)^2   E(100GeV)=mc^2
  double rm = 4*Mwimp*Mt/pow(Mwimp+Mt,2); 

  //  double Term1 = sqrt(pi)*v0/(4*vear) * ( TMath::Erf((vmin+vear)/v0) - TMath::Erf((vmin-vear)/v0) );
  double Term1 = c1 * exp(-c2 * ER / E0 / rm);
  double Term2 = exp( -pow(vesc/v0,2) );
  
  //  double dRdE = k0/kf * R0/E0/rm * ( Term1 - Term2 ); //* pow(F(Q,A,Mt),2);
  double drate =  R0/E0/rm * exp(-ER/E0/rm) *1e-6;

  return drate;
}

// reduce mass
double Mr(double M1, double M2){
  return M1*M2/(M1+M2);
}

//form factor  Helm form factor (CDMS thesis)           

double F(double ER, double A, double mT){ // ER -> recoil E , mT -> target mass (keV)
  int option=0;
  
  double Q = sqrt(2*mT*ER)/c; // keV s/m                                              

  double aT = 0.52                    *pow(10,-15); // m                              
  double cT = (1.23*pow(A,1/3.) - 0.6)*pow(10,-15); // m                              
  double sT = 0.9                     *pow(10,-15); // m                              

  double rT = sqrt( pow(cT,2)+7/3.*pow(pi*aT,2)-5*pow(sT,2) ); // [m] effective nuclear radius  

  double hbar = 6.582e-19; // keVs                                                    
  double xq = Q*rT/hbar; // dimentionless                                             
  
  double formfac = 3*(sin(xq) - xq*cos(xq))/pow(xq,3) * exp( -pow(Q*sT/hbar,2)/2. );       
  if(option==1) formfac = 3*sqrt(pi/(2*xq))*TMath::BesselJ1(xq)/xq * exp( -pow(Q*sT/hbar,2)/2. );  // use TMath BesselJ1 if you want 

 return formfac;
}

//-----------------                                                                   
// test form factor                                                                   
//-----------------                                                                   
double FormFac(double *x, double *par){
  double ER1 = x[0];
  double A1 = par[0];
  double mT1 = par[1];
  return pow(F(ER1,A1,mT1),2);
}

// for plot  -> dRdE * FormFactor

double dForm(double *x, double *par){
  double er = x[0];
  double A  = par[0];
  double sigmawn  = par[1];
  double mT = par[2];
  
  double er2=er * 1e-6;
  double drF = dRdE(er2,A,sigmawn) * pow(F(er,A,mT),2);
  //  double drF = dRdE(er2,A,sigmawn);
  return drF;
}

void drdeF(){
  TF1 *f1 = new TF1("f1", dRdE2, 0, 120, 2);
  //f1->SetParameters(40,1e-6); 
  f1->SetParameters(131,1e-6);     // A ->Xe , sigma(Wimp-n) =1e-6
  //  f1->SetParameters(130,1e-6); 

  TF1 *f2 = new TF1("f2", dForm, 0, 120, 3);
  f2->SetParameters(40,1e-6,40*0.932e6);    // A->Ar(40), sigma , mT =A*0.932(keV)
  //  f2->SetParameters(131,1e-6,131*0.932e6);
  TF1 *f3 = new TF1("f3", dForm, 0, 120, 3);
  //  f2->SetParameters(40,1e-6,40*0.932e6);  
  f3->SetParameters(131,1e-6,131*0.932e6);  
  TF1 *f4 = new TF1("f4", dForm, 0, 120, 3);
  //  f2->SetParameters(40,1e-6,40*0.932e6);
  f4->SetParameters(72.6,1e-6,72.6*0.932e6); // A->Ge(72.6)


  TCanvas *can = new TCanvas("can", "can");
  //  can->SetLogx();
  can->SetLogy();
  f1->SetTitle("dRdE");
  f1->SetLineColor(2);
  f1->Draw("");

  f2->SetLineColor(6);
  //f2->Draw("");
  f2->Draw("same");
  //f1->Draw("");
  f3->SetLineColor(1);
  f3->Draw("same");
  //  f2->Draw("same");
  //  f3->Draw("same");
  f4->SetLineColor(4);
  f4->Draw("same");

  f1->GetXaxis()->SetTitle("Recoil Energy [keV]");
  f1->GetYaxis()->SetTitle("Count [kev/kg/day]");
}

// plot for FormFactor by junji sanx
void drForm(){
  TF1 *f1 = new TF1("f1", FormFac, 0.00001, 500, 3);
  f1->SetParameters(40, 37.3e6, 0);
  TF1 *f2 = new TF1("f1", FormFac, 0.00001, 500, 3);
  f2->SetParameters(131, 122.7e6, 0);

  TF1 *f3 = new TF1("f3", FormFac, 0.00001, 500, 3);
  f3->SetParameters(40, 37.3e6, 1);
  TF1 *f4 = new TF1("f4", FormFac, 0.00001, 500, 3);
  f4->SetParameters(131, 122.7e6, 1);

  f1->SetNpx(10000);
  f1->SetMinimum(1e-5);
  TCanvas *c1 = new TCanvas("c1", "c1");
  gPad->SetGrid();
  gPad->SetLogy();
  f1->SetLineColor(2);
  f2->SetLineColor(1);
  f3->SetLineColor(2);
  f4->SetLineColor(1);
  f3->SetLineStyle(2);
  f4->SetLineStyle(2);
  f1->Draw();
  f2->Draw("same");
  f3->Draw("same");
  f4->Draw("same");

  f1->SetNpx(10000);
  f1->SetMinimum(1e-5);
  TCanvas *c1 = new TCanvas("c1", "c1");
  gPad->SetGrid();
  gPad->SetLogy();
  f1->SetLineColor(2);
  f2->SetLineColor(1);
  f3->SetLineColor(2);
  f4->SetLineColor(1);
  f3->SetLineStyle(2);
  f4->SetLineStyle(2);
  f1->Draw();
  f2->Draw("same");
  f3->Draw("same");
  f4->Draw("same");

  f1->GetXaxis()->SetTitle("Recoil Energy (keV)");
  f1->GetYaxis()->SetTitle("Squared Form Factor");
  f1->SetTitle("");

  TLegend *leg = new TLegend(0.7, 0.7, 1.0, 1.0);
  leg->AddEntry(f1, "Ar", "l");
  leg->AddEntry(f2, "Xe", "l");
  leg->AddEntry(f3, "Ar: Bessel", "l");
  leg->AddEntry(f4, "Xe: Bessel", "l");
  leg->Draw();
}
